9. Elastic Deformation#

Elastic deformation is computed using the ContactMechanics library.

Under the assumption that the surfaces are linear-elastic and isotropic, the deformation is computed using a boundary element method. Here, we are just interested in the displacement of the surface and not of the deformation within the material. This renders the 3D problem into a 2D problem that can be efficiently solved using fundamental solutions (so-called Green’s functions).

Theory#

Green’s Function Approach#

Assuming that the elastic deformation of the surfaces responds on a shorter timescale than the fluid-dynamic system, we demand the stresses \(\sigma\) in the material to be in equilibrium (\(\frac{\partial \underline{\sigma}}{\partial t} = 0\)), leading to the Navier-Lame equation.

\[ \nabla \cdot \underline{\sigma} = 0 \]

After putting in Hooke’s law relating stress to strain or displacement \(u\), respectively, we obtain a non-linear PDE in three dimensions. With some tricks, this system can be analytically solved for specific cases, see e.g. Johnson, Cambridge University Press, 2012. Contact Mechanics.

\[ \mu \nabla^2 \mathbf{u} + (\lambda + \mu) \nabla (\nabla \cdot \mathbf{u}) = 0 \]

To obtain the fundamental solution (Green’s function) for our approach, we solve the problem of a normal point force \(F_z\) acting on a infinite surface in \(x\)- and \(y\)-direction (infinite half-space) without any tangential or body forces. With this solution we can then calculate the normal displacement:

This is the essence of this approach: Solution of the displacement \(u_z(x,y)\) is obtained by convolution of the Green’s function \(G(x,y)\) with the pressure field \(p(x,y)\) over the domain \(A_G\). Note that \(G(x,y)\) only needs to be computed once.

\[\begin{split} \begin{aligned} u_z(x, y) & =\int_{A_G} G\left(x-x^{\prime}, y-y^{\prime}\right) p\left(x^{\prime}, y^{\prime}\right) d x^{\prime} d y^{\prime} \\ & =G(x, y) * p(x, y) \end{aligned} \end{split}\]

The Green’s function is different for periodic and non-periodic boundary conditions. See the following sections for details on

periodic solution
non-periodic solution
semi-periodic solution (periodicity in one direction, non-periodic in the other)

Periodic Solution#

For periodic boundary conditions, the system can be solved in Fourier space. The resulting Fourier space Green’s function takes the form:

\[ \tilde{G}_{\text{per}}(\vec{q}) = \frac{2}{E_{eff} \vert \vec{q} \vert} \]

with the effective Young’s modulus \(E_{eff}\) and the wavevector \(\vec{q}\). Note that we assume infinite thickness of the elastic half space. We set \(\tilde{G}(q_0) = 0\), resulting in the mean displacement being zero.

Non-Periodic Solution#

The non-periodic, real-space Green’s function takes the form

\[ G_{\text{nonper}}(x,y) = \frac{1}{\pi E_{eff}} \cdot \frac{1}{\sqrt{x^2 + y^2}} \]

with a singular point at the origin. To manage this singularity, the problem can be regularized by discretizing the domain into rectangles and to assume constant pressure over these rectangles.

By that, we obtain the real space solution for displacement of a general point \((x, y)\) due to a uniform pressure acting on a rectangular area \(2a \cdot 2b\) with center in the origin

\[ G_{\text{nonper}}(x, y) = \frac{1 - \nu^2}{\pi E} \int_{-a}^{a} \, \int _{-b}^{b} \frac{1}{\sqrt{(x - \xi)^2 + (y - \eta)^2}} \, d\xi\, d\eta \]

with the analytic solution

\[\begin{split} \begin{aligned} \frac{\pi E}{1-\nu^2} G_{\text{nonper}}(x, y) &= (x+a) \ln \left[\frac{(y+b)+\sqrt{(y+b)^2+(x+a)^2}}{(y-b)+\sqrt{(y-b)^2+(x+a)^2}}\right] \\ &+(y+b) \ln \left[\frac{(x+a)+\sqrt{(y+b)^2+(x+a)^2}}{(x-a)+\sqrt{(y+b)^2+(x-a)^2}}\right] \\ &+(x-a) \ln \left[\frac{(y-b)+\sqrt{(y-b)^2+(x-a)^2}}{(y+b)+\sqrt{(y+b)^2+(x-a)^2}}\right] \\ &+(y-b) \ln \left[\frac{(x-a)+\sqrt{(y-b)^2+(x-a)^2}}{(x+a)+\sqrt{(y-b)^2+(x+a)^2}}\right] \end{aligned} \end{split}\]

Semi-Periodic Solution#

This is a special case that might be needed, e.g. when simulating a section of a needle bearing with a 2D roughness profile, where we assume periodicity along the axis but non-periodic boundary conditions in the sliding direction.

Here, we numerically approximate the Green’s function by using the non-periodic solution and adding up periodic images in the direction where we want periodicity. For example, for periodicity in \(x\)-direction and non-periodicity in \(y\)-direction:

\[ G_{\text{semi}}(x,y) = \sum_{i=-n}^{n} G_{\text{nonper}}(x + i \cdot L_{x}, y) \]

with the length of the domain in \(x\)-direction \(L_x\) and the number of summed periodic images in both positive and negative direction \(n\). For our purposes, we assume \(n = 10\) providing sufficient accuracy.

Fourier Transform Trick#

Since the convolution operation scales with \(N^2\) in real space (\(N\) being the number of grid points), we exploit the fact that in Fourier space, convolution is just a “simple” matrix multiplication. As transforming to Fourier space and back-transforming to real space can be done efficiently using Discrete Fourier Transform (DFT), this workflow scales only with \(N \cdot \text{log}(N)\) Stanley, Kato, J. Tribol. 119, 481 (1997).

Note that for the non-periodic and semi-periodic case, the technique of decoupling periodic images is required to use this Fourier Transform Trick. For that, we assume that the pressure outside the domain is zero. Refer e.g. to Hockney, Methods Comput. Phys. 9, 135 (1970) and Pastewka, Robbins, Appl. Phys. Lett. 108, 221601 (2016).

Usage#

YAML#

If you want to include elastic deformation in your model, you need to introduce the elastic field under properties. Therein, you can specify the material parameters

Young’s modulus \(E\) in Pa, defaults to 210 GPa
Poisson’s ratio \(v\) in - , defaults to 0.30

as well as solver-related underrelaxation parameter that dampens the elastic deformation response. If you experience issues after enabling elastic deformation, try making this value smaller.

alpha_underrelax, defaults to 1e-03

In the special case of semi-periodic boundary conditions, you can set the number of periodic images (evaluated in both directions, i.e. n_images=10 translates to 21 images being summed up to obtain the Green’s function).

n_images, defaults to 10

properties:
    elastic:
        E: 210e09
        v: 0.3
        alpha_underrelax: 1e-03
        n_images: 10

Implementation Notes#

Periodicity is automatically derived from the density boundary conditions. If density is set periodic, elastic deformation will be calculated using the periodic Green’s function.
For non-periodic deformation, we define a reference pressure and reference deformation at the left side of the domain. This value is in both cases subtracted from all field values. Thereby, deformation is fixed to zero at the left boundary (This is a question of whether your given geometry should only define the initial state or if you want specific points to be fixed in the deformed state). See the Parabolic Slider, Non-Periodic example for clarification. You may want to change this fixation within the Topography.update() function within GaPFlow/topography.py.
Note that if you specify a non-periodic 1D problem in \(x\)-direction with \(N_y\)=1 and periodicity in \(y\)-direction, the solver will assume line contact for the calculation of the elastic deformation. That means that the non-periodic Green’s function is applied and \(L_y\) is set to the unit length of 1 m (This is relevant because the deformation is determined by the applied forces, and calculating the force requires specifying an area). As a consequence, this prevents accidental dependence of the result on \(L_y\).

Examples#

Journal Slider, Periodic#

This example shall highlight the deformation with periodic boundary conditions. Therefore we analyze a 1D journal slider and assume periodicity. Similar to the previous example, we assume a rather soft material with \(E = 5\) GPa to obtain significant deformation without requiring pressures in the GPa range.

sim = """
options:
    output: data/journal
    write_freq: 1000
    silent: False
grid:
    dx: 1.e-5
    dy: 1.
    Nx: 100
    Ny: 1
    xE: ['P', 'P', 'P']
    xW: ['P', 'P', 'P']
    yS: ['P', 'P', 'P']
    yN: ['P', 'P', 'P']
geometry:
    type: journal
    CR: 1.e-2
    eps: 0.7
    U: 0.1
    V: 0.
numerics:
    CFL: 0.25
    adaptive: 1
    tol: 1e-8
    dt: 1e-10
    max_it: 40_000
properties:
    shear: 0.0794
    bulk: 0.
    EOS: DH
    P0: 101325.
    rho0: 877.7007
    T0: 323.15
    C1: 3.5e10
    C2: 1.23
    elastic:
        E: 5e09
        v: 0.3
        alpha_underrelax: 1e-04
"""

Now we initialize the problem

from GaPFlow.problem import Problem

myProblem = Problem.from_string(sim)

*************************************************************
*                       PROBLEM SETUP                       *
*************************************************************
- options:
  - output                   : data/journal
  - write_freq               : 1000
  - use_tstamp               : True
  - silent                   : False
- grid:
  - Nx                       : 100
  - dx                       : 1e-05
  - Lx                       : 0.001
  - Ny                       : 1
  - dy                       : 1.0
  - Ly                       : 1.0
  - dim                      : 1
  - bc_xE_P                  : [True, True, True]
  - bc_xE_D                  : [False, False, False]
  - bc_xE_N                  : [False, False, False]
  - bc_xW_P                  : [True, True, True]
  - bc_xW_D                  : [False, False, False]
  - bc_xW_N                  : [False, False, False]
  - bc_yS_P                  : [True, True, True]
  - bc_yS_D                  : [False, False, False]
  - bc_yS_N                  : [False, False, False]
  - bc_yN_P                  : [True, True, True]
  - bc_yN_D                  : [False, False, False]
  - bc_yN_N                  : [False, False, False]
- geometry:
  - U                        : 0.1
  - V                        : 0.0
  - type                     : journal
  - flip                     : False
  - CR                       : 0.01
  - eps                      : 0.7
- numerics:
  - tol                      : 1e-08
  - max_it                   : 40000
  - dt                       : 1e-10
  - adaptive                 : True
  - CFL                      : 0.25
  - MC_order                 : 1
- properties:
  - shear                    : 0.0794
  - bulk                     : 0.0
  - EOS                      : DH
  - rho0                     : 877.7007
  - P0                       : 101325.0
  - C1                       : 35000000000.0
  - C2                       : 1.23
  - elastic:
    - enabled                : True
    - E                      : 5000000000.0
    - v                      : 0.3
    - alpha_underrelax       : 0.0001
    - n_images               : 10
- gp:
- db:
- md:
*************************************************************
*                  PROBLEM SETUP COMPLETED                  *
*************************************************************
                                                             
     Writing output into: data/2026-02-26_113452_journal

Let’s look at the initial geometry:

%matplotlib inline
myProblem.plot_topo()

../_images/5a03b33097eabe0f64029a965da00443e4e75042634a9a7fcd1475550a5536f4.png

Run the simulation and observe the result:

myProblem.run()

myProblem.plot_topo(show_defo=True, show_pressure=True)

-------------------------------------------------------------
Step   Timestep   Time       CFL        Residual  
-------------------------------------------------------------
    1.8984e-10 0.0000e+00 2.5000e-01 1.0000e+00
 1.8982e-10 1.8982e-07 2.5000e-01 3.4048e-04
 1.8981e-10 3.7964e-07 2.5000e-01 2.9804e-05
 1.8981e-10 5.6945e-07 2.5000e-01 7.9710e-06
 1.8981e-10 7.5927e-07 2.5000e-01 3.5214e-06
 1.8981e-10 9.4908e-07 2.5000e-01 1.5100e-05
 1.8981e-10 1.1389e-06 2.5000e-01 2.6030e-05
 1.8981e-10 1.3287e-06 2.5000e-01 3.5601e-05
 1.8981e-10 1.5185e-06 2.5000e-01 4.3432e-05
 1.8981e-10 1.7083e-06 2.5000e-01 4.9382e-05
1.8981e-10 1.8981e-06 2.5000e-01 5.3484e-05
1.8981e-10 2.0880e-06 2.5000e-01 5.5882e-05
1.8981e-10 2.2778e-06 2.5000e-01 5.6785e-05
1.8981e-10 2.4676e-06 2.5000e-01 5.6431e-05
1.8981e-10 2.6574e-06 2.5000e-01 5.5059e-05
1.8981e-10 2.8472e-06 2.5000e-01 5.2893e-05
1.8981e-10 3.0370e-06 2.5000e-01 5.0137e-05
1.8982e-10 3.2268e-06 2.5000e-01 4.6966e-05
1.8982e-10 3.4167e-06 2.5000e-01 4.3530e-05
1.8982e-10 3.6065e-06 2.5000e-01 3.9951e-05
1.8982e-10 3.7963e-06 2.5000e-01 3.6331e-05
1.8982e-10 3.9861e-06 2.5000e-01 3.2747e-05
1.8982e-10 4.1759e-06 2.5000e-01 2.9261e-05
1.8982e-10 4.3657e-06 2.5000e-01 2.5920e-05
1.8982e-10 4.5556e-06 2.5000e-01 2.2757e-05
1.8982e-10 4.7454e-06 2.5000e-01 1.9794e-05
1.8982e-10 4.9352e-06 2.5000e-01 1.7047e-05
1.8982e-10 5.1250e-06 2.5000e-01 1.4523e-05
1.8982e-10 5.3148e-06 2.5000e-01 1.2222e-05
1.8982e-10 5.5046e-06 2.5000e-01 1.0144e-05
1.8982e-10 5.6945e-06 2.5000e-01 8.2807e-06
1.8982e-10 5.8843e-06 2.5000e-01 6.6248e-06
1.8982e-10 6.0741e-06 2.5000e-01 5.1654e-06
1.8982e-10 6.2639e-06 2.5000e-01 3.8906e-06
1.8982e-10 6.4537e-06 2.5000e-01 2.7878e-06
1.8982e-10 6.6436e-06 2.5000e-01 1.8435e-06
1.8982e-10 6.8334e-06 2.5000e-01 1.0446e-06
1.8982e-10 7.0232e-06 2.5000e-01 3.7756e-07
1.8982e-10 7.1470e-06 2.5000e-01 7.4811e-09
=================================
Total walltime   :  0:02:41
(232.95 steps/s)
=================================

../_images/608c4b1abe0e91363637553762259a570d0c3ccf512895aa9aed53c413c66334.png

We can also generate an animation of the solution process (make sure you have the ipympl package installed):

myProblem.animate()

Note that, by using the Dowson-Higginson equation of state, we allow negative pressures in this example. For more physical behaviour, we should use an equation of state that considers cavitation and therefore does not output negative pressures. For that, see the non-periodic example below that uses Bayada-Chupin equation of state.

Parabolic Slider, Non-Periodic#

In this example, we want to illustrate how a parabolic slider deforms.

For this notebook, we will define this problem using a string variable as input. The structure is similar to the .yaml file. This is a 1D problem with non-periodic boundary conditions, therefore the elastic deformation will be computed using the non-periodic Green’s function. Note that we use a rather soft material with \(E = 5\) GPa to obtain significant deformation without requiring pressures in the GPa range.

sim = """
options:
    output: data/pslider1d_elastic
    write_freq: 1000
    use_tstamp: True
grid:
    Lx: 0.0762
    Ly: 1.
    Nx: 100
    Ny: 1
    xE: ['D', 'N', 'N']
    xW: ['D', 'N', 'N']
    yS: ['P', 'P', 'P']
    yN: ['P', 'P', 'P']
    xE_D: 850. # 101325
    xW_D: 850. # 101325
geometry:
    type: parabolic
    hmin: 0.5e-5
    hmax: 5.0e-4
    U: 10.0
    V: 0.
numerics:
    adaptive: 1
    CFL: 0.45
    tol: 1e-8
    dt: 1.e-10
    max_it: 100_000
properties:
    EOS: Bayada
    rho0: 850.
    shear: 0.039
    bulk: 0.
    cl: 1600.
    cv: 352.
    elastic:
        E: 5e09
        v: 0.3
        alpha_underrelax: 1e-03
    piezo:
        name: Dukler
        shearv: 3.9e-5
        rhol: 850.
        rhov: 0.019
"""